Instability of synchronized motion in nonlocally coupled neural oscillators 
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We study nonlocally coupled Hodgkin-Huxley equations with excitatory and inhibitory synaptic 
coupling. We investigate the linear stability of the synchronized solution, and find numerically 
various nonuniform oscillatory states such as chimera states, wavy states, clustering states, and 
. . , spatiotemporal chaos as a result of the instability. 
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Synchronization of neural oscillators is considered to play an important role in various functions such as visual 
OO • information processing, sleeping, and memory in the brain [^^|3|- ^ ne synchronization of neural oscillators has 
been theoretically and numerically studied by many authors la 13 • The conditions of the synchronization in 
the neural oscillators were intensively studied, but it is not well known what happens after the synchronized motion 
becomes unstable. Synchronization is not always desirable; several neurological diseases such as Parkinson's disease 
and epilepsy are caused by synchronized firing of neural oscillators. The control and prediction of synchronized motion 
are related to therapies for the diseases. A therapy for Parkinson's disease is external electric stimulation at high 
frequencies, called deep brain stimulation 0. Time series analyses of electroencephalograms were studied to predict 
an epilepsy seizure beforehand |10| . Some drugs for neurological diseases are interpreted to modify synaptic currents, 
change the interaction among neural oscillators, and make anomalously synchronized states change into asynchronous 
states. 

The stability of the synchronized solution depends on the details of the model equations of neural oscillators. 
t-H . However, there is some general tendency that excitatory coupling induces synchrony and inhibitory coupling leads to 
antisynchrony (in which the phase difference between two oscillators is ir) , if the response time of synaptic couplings is 
rather small. However, if the response time of synaptic couplings is slower, the situation is reversed. The synchronous 
state is stable for inhibitory coupling, and the synchronous state becomes unstable for excitatory coupling. The 
response time or the delay time is an important parameter for the stability of the synchronized motion. The stability 
of the synchronous state has been studied mainly in coupled two-neuron systems or in uniformly or randomly coupled 
\Q ' systems of many neurons. On the other hand, Kuramoto and collaborators found various complicated solutions 
in nonlocally coupled complex Ginzburg-Landau equations and nonlocally coupled phase oscillators 0, [^. We 
studied nonlocally coupled noisy integrate-and-fire models using the Fokkcr-Planck equation and found traveling wave 
solutions [l3| . In this paper, we study the instability of the synchronized solution in nonlocally coupled Hodgkin- 
i-^p Huxley equations by changing the response time of the synaptic coupling. 
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II. NONLOCALLY COUPLED HODGKIN-HUXLEY EQUATION 
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The Hodgkin-Huxley equation is a fundamental model equation for the firing of neurons. The equation was originally 
proposed for the giant axon in a squid [l4|. Generalized Hodgkin-Huxley type equations are used for various neurons. 
These are coupled differential equations for the membrane potential V and the auxiliary variables m, n, h, which are 
related to the conductance of ion currents through the Na + and K + channels and a leak current. The Hodgkin-Huxley 
equation is written as 
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where C m = 1, jNa = 120, <?k = 36, <?l = 0.3, i?Na = 50, Ek = —77 and — —54.4, I denotes the external current 
stimulus, and I s is the synaptic current. The unit of electric potential is the millivolt and the unit of time is the 
millisecond. We propose a time evolution equation for the synaptic current I s in a one-dimensional system as 

dl f°° 

T ^ = J 9(x-x')F(V(x'))dx'-I s , (2) 

where g{x) represents a coupling function of the nonlocal synaptic interaction, r is the response time, and F(V) is 
the response function of the synaptic current to the membrane potential. In this paper, we assume a function 

F(V) = 7 (V - V ) for V > V , 

= for V < V , (3) 

where the parameters Vq = —50 and 7 = 0.01 are used in this paper. As a coupling function g(x), we use an 
exponential function g(x) — g\a/2 exp(— ot\x\) or a sum of two exponential functions g(x) = giai/2exp(— ai|a;|) + 
giotil'Z exp(— o^M) for the sake of simplicity. For synchronized solutions, V(x, t), m(x, t), h(x, t) and n(x,t) do not 
depend on x; then the uniform solutions Vo(t), mo(t), ho(t) and no(t) obey the same equation (1), but the synaptic 
current I s o obeys 

r d -^=g F{V )-I s , (4) 

where g = g(x)dx. 

III. INSTABILITY OF SYNCHRONIZED MOTION IN EXCITATORY AND INHIBITORY SYSTEMS 

First, we investigate the linear stability of the uniform solution for the nonlocally coupled Hodgkin-Huxley equation. 
Small perturbations with wave number k around the uniform solution: SV = V k exp(ikx) , 5m — m k exp(ikx),5h = 
h k exp(ikx),5n = n k exp(ikx) and SI S = I sk exp(ikx) obey 

Cm ~dT = ffNa { 3m o m fcfro(£Na - v o) + mlh k (E Na - V Q ) - m^oVfc} 
+ g K 4nln k (E K - V) + g K n%(-V k ) + g^{-V k ) + I sk , 
dm k da m ( da m df3 m \ 
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T ^r = 9 k -gTrV k -i sk , (5) 

where g k is the Fourier transform of g(x). If the uniform solution Vo, m , h , n , and I s0 is time periodic, the small 
pertur bation increases or decay s periodically. We can calculate numerically the linear growth rate of the norm 
N = \JV k + m k + h? k + n k + I^ k . Numerical simulations are performed using the Runge-Kutta method. The system 
system size is L = 512 and the grid sizes At = 0.005 and Air = 1, and periodic boundary conditions are imposed. 
That is, 512 oscillators are set on a one-dimensional circle of length 512. Figure 1(a) displays the average linear 
growth exponent A = lim(l/T)ln[iV(T)/iV(0)] of the norm for large T as a function of k for an excitatory coupling 
g{x) = 1.8 cxp(— 0.03|a;|) and the uniform input Iq = 15. The four curves correspond to four response times r = 6, 7, 8 
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FIG. 1: (a) Linear growth exponent A for the perturbation with wavenumber k for the excitatory coupling g(x) 
1.8exp(— 0.03|a;|). (b) A for k = -n as a function of the response time r. 
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FIG. 2: Raster plots of firing neurons at r = 8.5 for (a) 2500 < t < 2550, (b) 3500 < t < 3550 and (c) 4950 < t < 5000. 



and 10. The synchronized state is stable at r = 6 and unstable at r = 7 for all wave numbers. Figure 1(b) displays the 
linear growth rate A for k = n. The instability occurs at t = 6.7. Figures 2(a), (b), and (c) display raster plots of firing 
neurons for some time intervals (a) 2500 < t < 2550, (b) 3500 < t < t < 3550, and (c) 4950 < t < 5000 at r = 8.5. 
In the raster plots, positions of neural oscillators satisfying the firing condition V(i,t) > are plotted with dots at 
every time interval 0.25 (ms). The initial condition is a uniform state with small random perturbation. The uniform 
state is unstable and the random perturbation enlarges in time as shown in Fig. 2(a). Discontinuities appear in the 
profile of membrane potential V(i), because the nonlocal coupling does not prohibit the discontinuity in contrast to 
the diffusion couplings. There appear mainly two discontinuous regions at i ~ 230 and i ~ 400 in Fig. 2(b). One 
discontinuous region survives at t = 5000 as seen in Fig. 2(c). The discontinuous region plays a role of a pacemaker 
and excitation pulses propagate toward both sides. Figure 3(a) displays a profile of the membrane potential V at 
r = 8.5. The membrane potential V is discontinuously distributed in 200 < i < 280. The time evolution of V is almost 
periodic, and the synaptic current I s changes in time almost periodically between 6 and 13. Figure 3(b) displays the 




FIG. 3: (a) Snapshot of the membrane potential V(i) at r = 8.5. (b) Frequency profile u)(i). 
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frequency profile u>(i) for the neural oscillators. The frequency of the oscillation in the discontinuous region is slightly 
faster than that in the continuous region, and it changes smoothly in space. This is similar to chimera states studied 
in Ref. [12,15]. Although they found the chimera states in conditions where the uniform state is stable, our chimera 
state has spontaneously appeared as a result of the instability of the uniform state. 
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FIG. 4: (a) Linear growth exponent A as a function of wave number k for the inhibitory coupling g(x) — — 1.8 exp(— 0.03js|). 
(b) Linear growth exponent A at k = it as a function of the response time r. 
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FIG. 5: (a) Raster plot for r = 4. (b) Time evolution of V(i) for i = 1 (solid curve) and i = 5 (dashed curve, (c) Snapshot of 
the membrane potential V(i) at t = 9985. 

Next, we consider an inhibitory system. The coupling function g(x) is assumed as g(x) — — 1.8 exp(— 0.03|cc|), and 
the uniform input is Iq — 15. Figure 4(a) displays the average linear growth rate A as a function of wavenumber k 
for r = 0.05,2,4, and 5. Figure 4(b) displays the average linear growth rate A for k = tt as a function of r. The 
uniformly synchronized state is unstable for r < 4.4. For the inhibitory coupling, the uniform state changes into a 
spatially modulated clustering state. Figure 5(a) displays a raster plot of firing neurons at r = 4. Figure 5(b) displays 
time evolutions of V(i) for i = 1 and 5. The time evolution is almost periodic, and it is clearly seen that the phase 
difference between the two oscillators i = 1 and 5 is about tt. Figure 5(c) displays a profile of the membrane potential 
V at t = 9985. The phase of each oscillator in each cluster is not uniform, but changes slowly in space. However, 
this phase profile is periodically recovered as seen in the raster plots in Fig. 5(a). The difference of the membrane 
potential V(i) between two oscillators inside the same cluster changes periodically but recovers the same value after 
one period of the oscillation. That is, the phase profile is random but it is frozen like a glassy state. 



IV. COMPETITION OF EXCITATORY AND INHIBITORY INTERACTIONS 



We can consider a more complicated system with excitatory and inhibitory interactions. We study the case of 
g(x) = — 1.8 exp(— 0.03|a;|) + 12 exp(— 0.12|a;|). The uniform input Iq — 15. The interaction is excitatory in the short 
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FIG. 6: (a) Linear growth exponent A for the perturbation with wave number k for the coupling g(x) = — 1.8 exp(— 0.03|x|) + 
12exp(— 0.12|a;|). (b) Raster plot to display firing neurons at r = 0.01. (c) Raster plot at r = 5. 
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FIG. 7: (a) Linear growth exponent A for the perturbation with wave number k for the inhibitory coupling g(x) = 
— 6exp(— 0.03|a;|) + 40 exp(— 0.12|a;|). (b) Raster plot to display firing neurons at r = 0.05. (c) Raster plot at r — 5. (d) 
Raster plot at r = 8. 



range and inhibitory in the long range. This type of interaction is often used in artificial neural network models. The 
uniform input is assumed to be Iq = 15. Figure 6(a) displays the average linear growth rate A as a function of wave 
number k for r = 0.01,5,6, and 7. The growth rate A becomes zero at k = k c = 0.139 for every r, where k c is a 
solution of g k = -1.8 • 0.03/(fc 2 + 0.03 2 ) + 12 • 0.12/(/c 2 + 0.12 2 ) = 40 = g , where g k is the Fourier transform of g(x). 
For r < 6.4, the growth rate A > for smaller k with k < k c , and A < for larger k with k > k c . On the other hand, 
in the case of r > 6.4, A < for k < k c , and A > for k > k c . There is a peak in the curve of A(fc) at a finite k, which 
corresponds to the most unstable mode. The existence of the characteristic wavelength is a result of the competition 
of the excitatory and inhibitory interactions. Figure 6(b) displays a raster plot of firing neurons at r = 0.01. A 
pacemaker region appears near i = 450. Traveling pulses are sent out regularly toward both sides from the pacemaker 
region. The spatial continuity is maintained for this coupling and therefore the frequency distribution is uniform 
in contrast to the chimera state. Figure 6(c) displays a raster plot of firing neurons at r = 5. The traveling wave 
state becomes unstable, spatial modulations grow up and a spatiotemporal chaos appears. Discontinuities sometimes 
appear in the raster plot as a result of the spatiotemporal chaos. 

Figure 7 shows a numerical result for a system with another coupling function g(x) = — 6 exp(— 0.1|a;|) + 
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40 exp(— 0.4|x|). In this coupling function, the critical wave number is k c — 0.464. The critical value of t for 
the instability is almost the same t c = 6.4 as the case shown in Fig. 7. The growth rate takes a maximum at a 
finite parameter k ~ 0.16 for r <~ 6.4. Figure 7(b) displays a raster plot of firing neurons at r = 0.05. Regular 
traveling waves state appear. There are three pacemakers. Because the ratio of the typical wave number for the case 
of Figs. 6 and 7 is 0.139/0.464 = 0.3, it is natural that the number of pacemakers is three times the case of Fig. 6. 
Figure 7(c) displays a raster plot of firing neurons at r = 5. Spatial modulations grow up and a spatiotemporal chaos 
appears. Discontinuities sometimes appear, and they induce phase slips. As a result, the frequency profile lo(x) is 
almost uniform, but randomly distributed around the uniform value. That is, the complete entrainment is broken by 
the phase slips. Figure 7(d) displays a raster plot of firing neurons at r = 8. Perturbations with small wavelength 
increase as seen in Fig. 7(d), because the linear instability occurs only for large wave number k > k c . 



V. INSTABILITY OF SYNCHRONIZED MOTION IN TWO-LAYER MODELS 

In neural systems, the characteristic of the interaction (that is, excitatory coupling or inhibitory coupling) is usually 
determined by the property of the presynaptic neuron. This is called Dale's principle. The interaction considered in 
the previous section is not so suitable in realistic neural systems. We can consider a two-layer model, where excitatory 
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FIG. 8: (a) Linear growth exponent A for the perturbation with wavenumber k for the two-layer model with coupling functions 
gi(x) = 5.4exp(— 0.12|rc|), 52 (z) = — 0.9 exp(— 0.03|x|), ; gz{x) = 3.6 exp(— 0.12|x|) and g±{x) = 0. (b) Maximum value of A 
(solid curve) and A^ (dashed curve) for the wavenumber k = 7v as a function of r. (c) Raster plot to display firing neurons in 
the first layer at r = 0.05. (d) Raster plot to display firing neurons in the first layer at r = 2. 

neurons in the first layer interact with inhibitory neurons in the second layer. The Hogdkin-Huxley equations for the 
two-layer model are written as 

^-JT = flN a m?/ii(£ Na -Vi) + g K ni(E K - Vi) + g L {E^ - Vi) + I Q1 + I sl + I s2 , 



dt 

Gm ~^ = • 9Nam 2 /l2 ( £ ' Na ~ + 9Knt(E K - V 2 ) + gL(E L - V 2 ) + I02 + ls3 + IsA, 



Tl ^T = J 9i{* - J)F(Vi{x'))dxf - I al , 
T2 ^df = P 92{x-x')F{V 2 {x'))dx' -I s2 , 
T3 ^df = I" 9a(x - x')F(Vi{x'))dx' - I a i, 
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FIG. 9: a) Linear growth exponent at r = 22,18 and 16 for the coupling functions gi(x) — 1.35 exp(— 0.03|»|), gi(x) = 
—3.6 exp(— 0.12|a;|), gz(x) = 3.6exp(— 0.12|:r|) and gi(x) = 0.(b) Raster plot of the neurons in the first layer at r = 22 



jr roc 

r 4 ^i = J g 4 (x - x')F(V 2 (x'))dx' - I sA , (6) 

where V\ and V% are membrane potentials of neurons in the first and second layers, and I s \, I S 2, I S 3, and I S 4 are, respec- 
tively, synaptic currents by the interaction inside the first layer, from the second layer to the first layer, from the first 
layer to the second layer, and inside the second layer. By Dale's principle, we assume that g\ > 0, <?2 < 0, 33 > and 
94 < 0. The interaction functions are assumed to be gi{x) = 5.4 exp(— 0.12|x|), 52(2^) = —0.9 exp(— 0.03|x|), 33(2;) = 
3.6 exp(— 0.12|x|) and g^ix) = 0. The response times n, T2, T3 and T4 are assumed to take the same value t for the sake 
of simplicity. These functions imply that the interaction among the neurons in the first layer is excitatory in the short 
range but the effective interaction becomes inhibitory in the long range via the inhibitory neurons in the second layer. 
The uniform inputs are assumed to be Iqi — 15 and /02 = 0, that is, the neurons in the second layer do not oscillate 
without the interaction with the first layer. Figure 8(a) displays the average linear growth rate A& as a function of k 
for r = 0.05, 1, 3 and 4. The behavior of the growth rate is more complicated, compared to the previous one-layer 
models. The linear growth rates are the largest at r = 1, that is, the synchronized state is most unstable for the 
intermediate response time. The wave number for the most unstable mode decreases as r is increased from 1 to 3, and 
becomes for r ~ 3.2. Figure 8(b) displays the maximum value of A and A w for the wavenumber k = tt as a function 
of r. The maximum growth rate increases with r in the small range of t, and then decreases for r > 1.2, and becomes 
for t > 3.2. The uniform state is stable for r > 3.2. Inhibitory interaction is dominant in this coupling function, 
because [gi(x] + g 2 {x)\dx is negative; therefore the larger response time makes the synchronized motion more 
stable. The linear growth rate for k — tt is positive for 0.61 < t < 2.1. However, the critical response time t c = 3.2 is 
rather smaller compared to the critical response time r = 6.4 of the one-layer model as shown in Figs. 6 and 7. This 
means that the two-layer system is favorable for the synchronized motion, which is consistent with previous studies 
@. Figure 8(c) displays a raster plot of firing neurons in the first layer at r = 0.05. Clustering domains appear, and 
the phase difference between neighboring clustering domains is nearly tt. Most neighboring neurons belongs to the 
same domain and discontinuities appear only between the neighboring domains, in contrast to the purely inhibitory 
system shown in Fig. 5. The sizes of the clustering domains are randomly distributed, and they depend on the initial 
conditions, but they are roughly determined by the wavelength I ~ 27r/0.06 ~ 105 of the most unstable Fourier mode 
in Fig. 8(a). Figure 8(d) displays a raster plot of firing neurons in the first layer at r = 2. A traveling wave state 
with two pacemakers appears. 

We have investigated another two-layer model with 31(2;) = 1.35 exp(—0.03|a;|), 32(2;) = —3.6 exp(— 0.12|a;|), gz(x) = 
3.6 exp(— 0.12|x|), and g±(x) = 0. For this coupling, the interaction among the neurons in the first layer is excitatory 
in the long range, but the effective interaction is inhibitory in the short range via the inhibitory neurons in the second 
layer. The excitatory coupling is dominant in this model, because J_ \jjx{x) + g2(x)]dx is positive. Figure 9(a) 
displays the average linear growth rate A& for t = 22, 18 and 16. At r = 16, the synchronized state is stable. At 
t = 18, the linear growth rate is weakly positive only for very small wave number. At r = 22, the linear growth rate 
is positive for all wave numbers. The weak instability for very small wave number occurs at r = 16.9 and the linear 
growth rates for all wave numbers become positive at r = 21.7. In any case, the critical response time is rather larger 
than the case of the one-layer model shown in Fig. 1. This is another example where the interaction via inhibitory 
interneurons makes the synchronization more stable. Figure 9(b) displays a raster plot of the firing neurons in the 
first layer at r = 22, where the short-wavelength instability occurs. The chimera state appears again also in this 
model as a result of the instability of the synchronized motion. 
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VI. SUMMARY 



We have studied the linear stability of the synchronized motion in the nonlocally coupled Hodgkin-Huxley equations. 
We have found various nonuniform oscillatory states such as chimera states, wavy states with several pacemakers, 
antisynchronous clustering domains, and spatiotemporal chaos as a result of the instability of synchronized motion. 
We did not show numerical results, but qualitatively the same results were obtained for a coupling function of a sum 
of two Gaussian functions g(x) = g\Oi\j ^Jlx exp[— (ai2p2] + §2012/ y/ft enp[—(oe2x) 2 ]. The nonlocal coupling is essential 
for the appearance of discontinuities in the chimera states, clustering states, and spatiotemporal chaos with phase 
slips. 

Some types of asynchronous states might be useful for information processing. It is important to make anomalous 
synchronization change into an asynchronous state in the case of neurological diseases. We hope that our numerical 
results of various asynchronous states might be relevant to these applications. 
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